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We develop a continuum theory to model low energy excitations of a generic four-band time 
reversal invariant electronic system with boundaries. We propose a variational energy functional 
for the wavefunctions which allows us derive natural boundary conditions valid for such systems. 
Our formulation is particularly suited to develop a continuum theory of the protected edge/surface 
excitations of topological insulators both in two and three dimensions. By a detailed comparison 
of our analytical formulation with tight binding calculations of ribbons of topological insulators 
modeled by the Bernevig-Hughes-Zhang (BHZ) hamiltonian, we show that the continuum theory 
with the natural boundary condition provides an appropriate description of the low energy physics. 
As a spin-off, we find that in a certain parameter regime, the gap that arises in topological insulator 
ribbons of finite width due to the hybridization of edges states from opposite edges, depends non- 
monotonically on the ribbon width and can nearly vanish at certain "magic widths" . 
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I. INTRODUCTION 

One of the physically observable phenomena in topo- 
logical insulators (TI) is the existence of the linearly 
dispersing gapless edge (in two dimension (2D) or sur- 
face (in three dimension (3D)) states which are topolog- 
ically protected against moderate electronic interactions 
or nonmagnetic disorder. 1-5 The fact that these conduct- 
ing edge states can host spin current without dissipation 
makes topological insulators (TI) a promising candidate 
in technological applications. 6,7 Understanding the na- 
ture of the edge states has been an aspect of interest in 
theoretical studies of topological insulators. 8-12 

Properties of edge states can be studied by construct- 
ing appropriate tight binding model Hamiltonians of TIs 
and examining their eigenstates for lattices with bound- 
aries. An alternative route is to construct a low energy 
continuum theory 3 ' 9 ' 13 that allows for analytical treat- 
ment that aids the development of field theoretic de- 
scription in the presence of interactions. 14 ' 15 Such ap- 
proaches have been gainfully employed earlier in studies 
of graphene. 16-20 In the analytic calculation, the edge 
states are obtained by subjecting appropriate bound- 
ary condition (BC) on the wavefunction. Here one 
usually 8-12 imposes the fixed boundary condition (also 
called as essential or Dirichlct boundary condition in the 
mathematical literature 21 ) where the wavefunction is as- 
sumed to be zero at the boundaries or at a fictitious layer 
of atoms just outside the boundaries. Such a choice of 
BC constraints the nature of the wavefunction in that the 
maximum weight of the edge state does not occur in the 
edge layers but in bulk layers that are near the edge layer. 
In the presence of interactions, the edge states and bulk 
states mix and the ensuing physics is determined cru- 
cially by this mixing. In a recent study 22 , it was shown 
that the Mott transition in topological insulator ribbons 
can occur in two different routes - the synchronous and 
asynchronous routes - depending on the nature of edge 



states. A continuum field theoretic analysis of such a 
phenomenon, therefore, requires a careful treatment of 
the edge states so that their profile correctly captures 
the mixing with the bulk states. 

With this motivation, in this paper, we develop a con- 
tinuum theory of time reversal invariant four-band model 
Hamiltonians that have been extensively used in the anal- 
ysis of topological insulators in two and three dimensions. 
We construct an energy functional of the wave functions; 
the wave function that renders this energy functional ex- 
tremum is shown the satisfy a stationary Schrbdinger 
equation that matches the four-band lattice theory at 
long wavelengths. As a key outcome of this approach, 
we derive a new boundary condition, the natural bound- 
ary condition. 21 This boundary condition is valid for any 
four-band time reversal invariant system in two and three 
dimensions. We use the BHZ model 3 that has been stud- 
ied earlier 8-11 to show that wthin a regime of parameters 
of this model, the natural boundary condition provides 
an excellent description of the edge states. In the pro- 
cess of this study, we show that the gap that arises from 
the hybridization of the edge states localized on the op- 
posite edges of a ribbon is a non-monotonic function of 
the ribbon width. This finding could potentially be use- 
ful in many applications such as design of thermoelectric 
devices etc. 23 ' 24 

In the following section (sec. II) we introduce a gen- 
eral four-band lattice Hamiltonian that is time reversal 
invariant. Sec. Ill contains the continuum theory of these 
systems, the formulation of a variational principle and 
derivation of the boundary conditions. A detailed com- 
parison of the numerical tight binding calculations and 
the analytical continuum theory is carried out in sec. IV 
using the BHZ model, 3 in its topological regime. The 
paper is concluded in sec. V which contains a discussion, 
significance and summary of the results. 
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II. FOUR-BAND TIME REVERSAL 
INVARIANT SYSTEMS 

Consider a Bravais lattice in two or three dimensions 
whose sites are labelled by /. Each lattice site has two 
orbitals (or "basis" sites such as A-B sites in the graphene 
lattice, sometimes also referred to as "flavours") labelled 
by a. The operator C\ aa creates an electron of spin a 
(quantized along some convenient axis) in the orbital a 
at site /. The Hamiltonian of the system is given by 

IS 

where S runs over lattice vectors, summation over re- 
peated orbital and spin indices is implied. The hopping 
matrix elements t ac ,^ a i {8) are such that the Hamilto- 
nian eqn. (f) is time reversal invariant. Hamiltonians 
discussed in the literature on topological insulators 4,5 are 
of this type. 

With the aim of developing a long wavelength contin- 
uum theory of such systems, we cast the Hamiltonian in 
the reciprocal space: 

H = ]T H ab {k)Ct a C kb (2) 
fceB 

where a (and b) is an index that represents the composite 
aa. Repeated a and b indices are summed over and k 
runs over B, the Brillouin zone which is a torus for 2D 
systems and a 3-torus in 3D systems. Following Refs. [1, 
25, and 26], we now write the matrix H{k) in a basis 
of sixteen 4x4 matrices, broken up into two groups T m 
(m = - 5) and A" {n = 1 - 10), i.e, 

5 10 

H{k) = d n {k)T n + £ e m {k)A m (3) 

n— m— 1 

where d n {k) and e n {k) are smooth functions of k. The 
matrices T and A are defined using r and er, the 2x2 
Pauli matrices associated with the orbital and spin de- 
grees of freedom, and 1, the 2x2 identity matrix. We 
have, r° = 1 <x> 1, 

r l,2,3,4,5 = { T * (g>t,T Z (g)t,T V (g)(T X ,T V ®(T y ,T y (g) CT*} 

(4) 

The ten elements A™ 1 , m = 1, . . . , 10 can be obtained 
from the commutators [r™, F n }/ {2i),n — 1, . . . , 5, n' > 
n. The grouping of these matrices into Ts and As is mo- 
tivated by the fact that under the action of the time 
reversal operator = — i(l <g) a y )K where K is the 
complex conjugation operator 27 , _1 r™@ = T n while 
1 A m © = -A m . From the fact that the Hamiltonian 
in eqn. (2) is time reversal invariant, and from the prop- 
erties of the r and A matrices just mentioned, we get 
from eqn. (3) that 1 

d„{-k) = d„{k) ^ 



Eqn. 2 along with eqns. 3 and 5 describes a general four- 
band Hamiltonian with time reversal symmetry. 

The systems of interest are those which possess a gap in 
their energy dispersion - two bands and separated from 
the other two by an energy gap - and the chemical po- 
tential lies in this gap. The nature of this insulating state 
(topological or trivial) is determined by the topological 
properties of the occupied bands and is characterized by 
the Z2 index. 1-3,28-30 While our formulation is applicable 
to any four-band system with time reversal symmetry, we 
shall focus on topological insulators which possess pro- 
tected edge/surface states. 

III. CONTINUUM THEORY, VARIATIONAL 
PRINCIPLE AND BOUNDARY CONDITIONS 

The continuum theory is developed by focusing on a re- 
gion of the Brillouin zone, i. e., specifically around the k- 
points which support low energy excitations. In the case 
of topological insulators with a bounding edge (or sur- 
face), the low energy excitations (i. e., excitations close 
to the chemical potential) usually occur near a time re- 
versal invariant momentum (TRIM). 4 TRIMs occur at 
the origin of the Brillouin zone, at the zone edges etc. 
In what follows, we shall develop the continuum theory 
focusing on the k = TRIM; generalization to any other 
TRIM of interest is straightforward. 

We discuss the continuum theory in the first quantized 
form. For our four-band model, the wave function is a 
four component vector function ip a {x) of the position x. 
We look to determine a Hamiltonian operator H that 
dictates the time evolution of ip a (x), i. e., 

iip a (x) = M ab ijj b {x) (6) 

where the dot represents time derivative and the repeated 
index b is summed over. We have set H = 1. To determine 
H, we expand the function d n {k) and e m {k) about k = 
up to quadratic order, which upon using eqn. (5) gives 

d n {k) = d n + kiS^jki 
6m{k) = 2A" l ki 

where the constants tensors S 1 ™ and vectors A™ are 
properties of the four-band system that characterize the 
dispersion near k = 0. We thus have 

H ab {k) » {d° n + S™hk 3 )T n ab + 2A™hA2 (8) 

where repeated n and m indices are summed over the 
ranges indicated in eqn. (3). H ;, can now be obtained as 
H a6 = H ab {ki -idi) where di = d/dxi, i. e., 

H ab = {d° n - S^d 3 )T n ab - 2iA?di AS. (9) 

which completes the discussion of eqn. (6). 

Consider now a region of space (in two or three dimen- 
sions) fl bounded by a boundary d£l (which may be an 



3 



edge or a surface). The stationary states at low energy 
are eigenstates of the continuum Hamiltonian H, i. e., 

U ab ^b{x) =E^ a (x) (10) 



where E is the energy eigenvalue, with appropriate 
boundary conditions for the four component wavefunc- 
tion ^ a (x) on d£l. 

To aid the determination of the boundary conditions, 
here we propose an energy functional associated with a 
four component wavefunction ^ a (a;): 



£[V(r), *(r)] = jf d d r{^:d^T n ab ^ h (d^S^d^) - i K^A^A + {d^lWK^b] E^ a ) (11) 

where E is a Lagrange multiplier that ensures that the wavefunction is normalized to unity. All repeated indices are 
summed over their appropriate ranges. We now show that the states that render this energy functional extremal are 
the stationary states of eqn. (10). Towards this end, upon varying \&* by S^f*, we get 



se 



d' 



= / d rf r(^:) f((d° - S^d )T n ab - 2iA?d i k a l - E6 ab ) *J + / d^r^l) k {S™T n ab d^ b + LA?A£* b ) 



<Kl 



(12) 



where we have used the divergence theorem and is the outward normal to the boundary d£l. 

I 



The extremality of £ necessitates that 

((d° - S^did^ - ^ATd.AZ - E6 ab ) * h = (13) 

in which is exactly the stationary Schrodinger equation 
of eqn. (10). Further on the boundary dfl, we have either 

SK = (14) 

which corresponds the fixed boundary condition where 
the values of the wavefunction \£ a is fixed (usually to 
zero), or 

n t (s^r: b d 3 * fc + iATKb^b) = (15) 

which is the natural boundary condition (note, again, 
that all the repeated indices are summed). We empha- 
size that this boundary condition is applicable to any 
time reversal invariant four-band system in two or three 
dimensions. In particular, the formulation is tailor made 
for the study of edge (surface) states of topological insu- 
lators. In the next section, we illustrate this framework 
by calculating (analytically) the edge states of a topolog- 
ical insulator described by the well known BHZ model 3 . 

IV. BHZ MODEL: COMPARISON OF 
CONTINUUM THEORY AND TIGHT BINDING 
RESULTS 

The BHZ model 3 describes 2D topological insulators 
realized in the HgTe/CdTe quantum wells. The tight 
binding version of the model is obtained by consider- 
ing four spin-orbit coupled orbitals- \s f), \p t) = 
I (p y + iPx) t), |s I), and |p |) = | (p y - ip x ) |) per site 



on a square lattice whose lattice spacing a is taken as 
unity. The model can be written as, 

Iacr ISaficr 

(16) 

where a,/3 = s,p and e Q denote the orbital energies. 
a =t, 4- an( l <5 is a nearest neighbour vector. The hopping 
matrix elements t a p{8<r) in the \sa), \pa) basis are given 
by, 

t(±xa) = ( ^ ± f^) , t(±ya) = ( % J%) 

\ ±f7 vi f pp J V^vi tpp J 

(17) 

where t ss , t sp , t pp are overlap integrals and a = +1 (—1) 
for spin \ (|). In the reciprocal space, as in eqn. (2), this 
Hamiltonian is described by matrices 

*<*>-(*?'*.<%) < 18 > 

where 

h(k\ — ( e s — %t s (cosk x + cosky) 2t sp (smk x — isinky) 
\ 2i sp (sin k x + i sin k v ) e p + 2t s (cos k x + cos k y ) 

(19) 

where we have set t ss = t pp = t s . Further defining eo 
such that e s = — (eo — 4i s ) and e p — (eo ~ 4t s ) we have 

H(k) = d 2 {k)T 2 + ei(fc)A 1 + e 2 (fc)A 2 (20) 
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in the form of eqn. (3), with r 2 = t z ® 1, A 1 = t x 
cr z ,A 2 = T y®\ and 

d 2 (k) = — e + 2i s (2 — (cosfc x + cosA: y )) 

ei(fc) = 2t sp smk x (21) 

e 2 (fc) = 2t sp smk y 

All other d(k)-s and e(fe)-s are zero. Note that here we 
have relabelled the m index in cqn. (3) for convenience. 
With this, focusing on the TRIM at k = 0, we get the 
continuum Hamiltonian operator as 

h = (-e - t s (d 2 x + a 2 )) r 2 - 2zt sp a c A 1 - 2it sp d y K 2 

(22) 

with d% = —e ,Sfj = t s 5ij,Al = t sp ,A 2 y = t sp ; all other 
d-s, S-s, As are zero. This Hamiltonian, upon setting 
t s = 1 has two scales, eo and t sp . When eo > 0, the 
system is in the topological phase; the remainder of the 
discussion considers only this case. The quantity t sp is a 
measure of the hybridization of the s and p orbitals and 
determines the "multi-componentness" of the wavefunc- 
tions. It must be noted that this model conserves the 
spin quantum number, i. e., the j" and J, spins decouple 
at the one particle level. 

In order to study the edge states of this model, we 
consider a geometry with — (—00,00) x (0, L), i. e., 
and infinitely long (along x-direction) ribbon of width L 
(terminated at y — and y = L, i. e, d£l — (y = 0) U (y = 
L)). When L — > 00, we get a half-space. 

Since the spins sectors decouple, we shall consider 
only the f-spin sector; the results of the J.-spin sector 
can be obtained by a time reversal operation. Exploit- 
ing the translational invariance along the x-direction, we 
write ty a (x,y) = e zkXl ^ a (y). For a given momentum k, 
the functions Q a (y) satisfy eqn. (10) with H given by 
eqn. (22): 



-e + t s (k 2 -0 2 ) 2t sp {k-d y ) 
2t sp (k + d y ) e -t(kl-dl) 



E 



- vj 
(23) 



Defining $ = 4> s + $ p , * = * s - * p , G{D) = G(d y ) = 
-e + t s (k 2 - a 2 ) and H(D) = H{d y ) = ~2t sp d y , we get 



G(D)* - H(D)V = (E- 2t sp k)§ 



(24) 



G{D)$ + H(D)<P = (E + 2t sp k)^ 
which leads to 

\G(D) - H(D)][G(D) + H(D)}$ = (E 2 - 4t 2 p fc 2 )$ 

(25) 

Assuming a trial solution Q(y) — e qv , we obtain the fol- 
lowing quartic equation for q, 

t 2 s (k 2 - q 2 ) 2 + 2(-e t s + 2t%)(k 2 - q 2 ) + (e 2 - E 2 ) = 

(26) 



which gives four solutions for q, gi 2 
which are given by, 



±<li, 93,4 = ±qn 



a 2 - k 2 
1i.ii — K 



(-e t s + 2t%) ± y/4t%(t% - e t s ) + t 2 E 2 
t 2 



(27) 

Therefore the general solution for $ and ^ are given by, 

*(y) = Aie qiv + A 2 e q2V + A 3 e q3V + A 4 e qiV (28) 
1 



Hv) 



E ~\~ ^tgpk 



{G(D)+H(D)}<S>(y) 



(29) 



where Ai-s are four constants. The complete solution for 
the wavefunction is given by, 



I s 



ikx 



1 



(30) 



The determination of the energy eigenvalue E and the 
constants .Ai-s requires the boundary conditions. The 
fixed boundary condition 9 eqn. (14) reads 



* S (0) = * P (0)=0 
* s (L) = * p (L)=0 



(31) 



while the natural boundary condition derived in eqn. (15) 
provides 

d*s 
dy 
d*„ 



t. 



ispfp = 



(32) 



t. 



dy 



tsp^s = 



on dfl i.e., at y — and y = L. 

In the remainder of the discussion t s is set to unity. 
It is useful to discuss the nature of the solution of q, be- 
fore proceeding to compare the analytical results with the 
numerical tight binding calculations. Note that the val- 
ues of q depends on the energy eigenvalue E (eqn. (27)). 
Since the k — corresponds to a TRIM, we expect pair 
(time reversal related) of topologically protected edge 
states at k = 0, and by the symmetry of the problem, 
we expect E = to be their energy eigenvalue. The 
values of q with E = are then determined by the pa- 
rameters eo and t sp , i. e., they are characterized by the 
same parameters that determine the "topology" of the 
system. Fig. 1 shows a plot of the qs as a function of 
the parameter e . We find that there are two regimes 
of eo, e < t 2 p , where there are four distinct real roots 
for qs and eo > t 2 p where gs are complex and appear 
in conjugate pairs. In the former regime, magnitudes of 
qi and q2 increase with increasing eo, while that of q 3 
and <74 decrease with increasing eo- In the latter regime, 
the real parts of q are unaffected, while their imaginary 
parts increase in magnitude. Clearly, the nature of the 
edge states for e < t 2 p is different from that for e > t 2 p . 
In the former case, the edge state wavefunction is non- 
oscillating and falls exponentially as the distance from 
the edge. In the latter case, the wave function also has 
an oscillatory part, and as we shall show later, this leads 
to quite interesting physics and possibilities. 
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FIG. 1. The dependence of the wavevectors q that determine 
the nature of the edge states on eo. Top: Real part of q. 
Bottom: Imaginary part of q. For eo < t 2 ap , the edge states 
are exponentially decaying, while for eo > t^ p , they have an 
oscillating character along with the exponential decay. 



A. Half Space 

Let us first consider a semi-infinite plane with its 
boundary at y = 0. Then the bounded solution for 4> 
is given by, 



<t>{y) = A 2 e 



II v 



-guv 



(33) 



The energy eigenvalues and the wave functions can 
be determined by imposing either the fixed boundary 
condition eqn. (31) or the natural boundary condition 
eqn. (32). After some simple algebra, it can be shown 
that for small k, 



E(k) = 2t sp k 



(34) 



a linear dispersion for the edge states, that is, remarkably, 
independent of which boundary condition is chosen. 

This value of E(k) can be now used to determine the 
constant coefficients ^2,4 and hence ^ s ,p- The profile 
of the wave function, of course, depends strongly on the 
boundary conditions. Fig. 2 shows a comparison of the 
results of the analytical formulation presented above with 
the two different boundary conditions and the wave func- 
tion obtained from numerical calculations with the full 
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FIG. 2. Comparison of the edge state wavefunctions obtained 
analytically by using the two different BCs with the corre- 
sponding result from the tight binding calculation. The wave- 
functions plotted are for k = 0; only the ^ s component is 
shown. 



tight binding model. Fig. 2(a) shows that for a value of 
t sp = 0.5, the wave function calculated from with the 
fixed boundary condition differs significantly from that 
of the tight binding results for points close to the edge 
(near y = 0) . The wavefunction with the natural bound- 
ary condition does not vanish at the boundary and has 
the expected exponential decay into the bulk. At large 
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FIG. 3. Energy dispersion of edge states of a BHZ ribbon of 
width L — 20. For the parameter values shown, the contin- 
uum theory with the natural boundary condition (eqn. (15)) 
reproduces the tight binding result more accurately. 

distances from the boundary the tight binding result for 
the wave function falls between the analytical results of 
the fixed and natural boundary conditions. This can be 
understood by noting that the fixed boundary condition 
kills the weight of the edge state near the boundary, and 
hence overestimates the weight of the wave function in 
the bulk. The effect is precisely the opposite with the 
natural boundary condition, where the weight in the bulk 
is underestimated compared to tight binding result. We 
now consider Fig. 2(b) which shows the comparison of 
the edge state wave function with t sp = 2, but still with 
eq < t 2 sp . In this case we see that the wavefunction deter- 
mined by the natural boundary condition not only closely 
reproduces the qualitative aspects of the tight binding 
solution, but is also in excellent quantitative agreement 
with it at large distance from the edge. Finally, in Fig. 2 
we show the comparison of the wave functions in the 
regime of parameters with e > t%- We see, again, that 
the analytical wave function obtained with the natural 
boundary condition more closely matches the results of 
tight binding calculation. 



B. Ribbons 

We now consider ribbons of finite width L. In this 
case, the edge states emanating from the edges at y = 
and y = L, overlap and hybridize rendering the system 
gapped (see Fig. 3). A stronger test of the validity of 
the continuum formulation and the correctness of the 
boundary condition can achieved by comparing the gap 
calculated using the analytical formulation with that ob- 
tained from the tight binding numerics. Fig. 4(a) shows 
the comparison of the calculated gaps as a function of the 
ribbon width L. In this regime of parameters the gap falls 
exponentially with the ribbon width as it is determined 
by the overlap matrix element of the two edge states em- 
anating from the opposite edges. Again, we see that in 
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FIG. 4. Comparison of the energy gap obtained analytically 
by using the two different BCs with the corresponding result 
from the tight binding calculation. 



this parameter regime, the tight binding gap lies between 
the fixed boundary condition result which is the largest, 
and the natural boundary condition value which is the 
smallest. This can be understood based on the result of 
the previous section. The weight of the edge state wave 
function in the bulk is overestimated by the use of the 
natural boundary condition and hence this gives rise to a 
larger gap owing to a larger overlap of the wavefunctions 
emanating from the opposite edges. For the same reason, 
the natural boundary condition underestimates the gap. 
For a larger value of t sp , the natural boundary condition 
is in better quantitative agreement with the tight binding 
results. This owes, again to the fact that wave function 
is better estimated by the natural boundary condition. 

Our final result pertains to the energy gap in ribbons 
with parameters in the regime eo > t 2 sv - Fig. 5 shows a 
plot of the gap as a function of the ribbon width in such a 
regime; we see that the gap is non-monotonic. Although 
the gap follows an exponential fall with increasing rib- 
bon width, there are "magic widths" at which the gap 
is very small; indeed our analytical results with the nat- 
ural boundary conditions does reproduce these features. 
The physics behind this phenomenon can be traced to 
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oscillatory in nature, the current brought about by one 
flavour can be reflected in another flavour channel. They 
may be contrasted with systems with a single component 
wave function such as in a simple "one component" tight 
binding model where the appropriate continuum bound- 
ary condition is that the vanishing of the wavefunction 
at the boundary. Topological insulators that are "deep" 
in their topological phase (such as a large eo and t sp ) are 
strongly "multi-component" in nature. For such systems 
the natural boundary condition is more appropriate. 



V. SUMMARY 



FIG. 5. Non-monotonic dependence of the energy gap with 
ribbon width L in the parameter regime eo > t%- The dashed 
line is a guide to the eye to show an overall exponential de- 
pendence on the ribbon width. 

the oscillatory nature of the edge state wave function in 
this parameter regime; for some particular widths of the 
ribbon, there is a "near destructive interference" of the 
wave functions emanating from the opposite edges that 
renders their overlap matrix element small resulting in 
a smaller gap. To the best of our knowledge, this is the 
first report of such physics in the BHZ model. We believe 
this is generic, and in fact, can find possible use in the 
design nano-scale devices with topological insulators. 

C. Discussion 

As is evident from our results, a continuum field the- 
ory with a natural boundary condition provides an ex- 
cellent description of systems with strong "component- 
mixing". In the case of the BHZ model, this will occur 
when t sp is large. Physically, in such cases a wave of "one 
flavour" can be reflected off a boundary as another fla- 
vor, and thus the wave functions do not have to vanish. 
This applies to the regime were the wave functions are 



In the paper, we have developed a continuum theory 
that is applicable to study four-band time reversal in- 
variant systems. We formulate a variational energy func- 
tional and show that the Schrodinger equation in the 
bulk is the Euler-Lagrange equation of this functional. 
This formulation allows us to obtain the natural bound- 
ary condition of the system. We have compared our an- 
alytical results with full tight binding calculation for the 
BHZ model for the half-space and finite ribbons. We 
show that in the interesting topological regime, the nat- 
ural boundary condition derived in this paper is more 
appropriate. We believe that our continuum formulation 
and boundary conditions will be useful in developing the- 
ory of devices and applications of topological insulators, 
and continuum theory modeling of experiments such as 
tunnelling from surface states. The non-monotonic de- 
pendence of the gap on the width of a BHZ ribbon is of 
particular interest; we believe such features are generic 
and can have numerous applications. 
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